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ABSTRACT 

In this series of papers we study the structure of the atomic to molecular transition in the giant 
atomic-molecular complexes that are the repositories of most molecular gas in galaxies, with the ulti- 
mate goal of attaining a better understanding of what determines galaxies' molecular content. Here 
we derive an approximate analytic solution for the structure of a photodissociation region (PDR) in 
a cloud of finite size that is bathed in an external dissociating radiation field. Our solution extends 
previous work, which with few exceptions has been restricted to a one-dimensional treatment of the 
radiation field. We show that our analytic results compare favorably to exact numerical calculations 
in the one-dimensional limit. However, our more general geometry provides a more realistic represen- 
tation than a semi-infinite slab for atomic-molecular complexes exposed to the interstellar radiation 
field, particularly in environments such as low-metallicity dwarf galaxies where the curvature and 
finite size of the atomic envelope cannot be neglected. For clouds that are at least 20% molecular we 
obtain analytic expressions for the molecular fraction in terms of properties of the gas and radiation 
field that are accurate to tens of percent, while for clouds of lower molecular content we obtain upper 
limits. As a side benefit, our analysis helps clarify when self-shielding is the dominant process in H2 
formation, and under what circumstances shielding by dust makes a significant contribution. 
Subject headings: ISM: clouds — ISM: molecules — molecular processes — radiative transfer 



1. INTRODUCTION 

In galaxies such as the Milky Way, where atomic 
and molecular phases of the interstellar medium coex- 
ist, molecular clouds re present the inner parts of atomic - 
molecular complexes (Elmcgrc en fc Elmegreenl Il987h . 
The bulk of the volume of the interstellar medium is filled 
with far-ultraviolct (FUV) photons capable of dissociat- 
ing hydrogen molecules, and this radiation field keeps the 
majority of the gas atomic. Gas that is predominantly 
molecular is found only in dense regions where a com- 
bination of shielding by dust grains and self-shielding 
by hydrogen molecules excludes the interstellar FUV 
field. These molecular regions are bounded by a pho- 
todissociation re gion (PDR) in which the gas is predom- 
inantly atomic (|Hollenbach fc Tielensl Il999t and refer- 
ences therein). 

To date most work on the structure of PDRs has 
been limited to one-dimensional geometries, includ- 
ing unidirectional or bidirectional beams of radia- 
tion impinging on semi-infinite slabs or purely ra- 
dial radiation fields strik i ng the surfaces of spheres 
(e.g. iFederman et al.lll979t Ivan Dishoeck fc Black! I198H 



Black fc van DishoeckJ [l987t ISternberdH9 88: Elmegree, 
19931: iDraine fc Bertoldil 119961: iHollenbach fc Tielen 



1999; iBrowning et al.ll2003l : lAlle n et al.1 12004) . For these 
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one-dimensional problems the literature contains both 
detailed numerical solutions and analytic approximations 
for the problem of radiative transfer and H2 formation- 
dissociation equilibrium. These approaches yield good 
results when the PDR is thin compared to the cloud 
as a whole, or for PDRs that are in close proximity to 
hot stars whose radiation and winds have compressed 
the PDR into a slab-like geometry. Many nearby well- 
studied PDRs, such as the Orion bar, fall into this latter 
category. However, the one-dimensional approximation 
is much less appropriate for giant clouds being dissoci- 
ated by the combined starlight of many distinct stars 
and star clusters, particularly when the atomic region 
constitutes a significant fraction of the total cloud vol- 
ume. The problem is especially severe in galaxies with 
low metallicities and interstellar pressures, where the pre- 
dominantly molecular parts of cloud complexes generally 
constitute a small part of the total mass and volume (e.g. 
iBlitz fc Rosolowskvl 12006). In this case one cannot ne- 
glect either the curvature of the PDR or the finite size of 
the molecular region, and a higher-dimensional approach 
is preferable. Previous work in one dimension is there- 
fore of limited use for the problem on which we focus: 
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determining the atomic and molecular content of galax- 
ies on large scales, under the combined effects of all the 
sources of dissociating radiation in that galaxy. 

Considerably less work has focused on higher- 
dimensional geometries, since these require a treatment 
of the angular dependence of the radiation field and its 
variation with position inside a cloud. As a result all 
treatments of two- or three-dimensional radiation fields 
to date are purely numerical, 
consider spherical clouds and 
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allow arbitrary geometries, but their method applies only 
in translucent clouds, and involves an approximate nu- 
merical integr a tion o f the transf e r equ ation. Similarly, 
iLiszt & Lucas] (|2000D and iLisztl (|2002D present models 
for PDRs in spherical clouds involving angular integra- 
tion over the radiation field in radial bins, coupled with 
a relaxation method to determine the H2 abundance at 
each radius. Neither of these approaches yield a simple 
analytic estimate of the size of the PDR or the molecular 
region, nor do they provide any insight into the dimen- 
sionless numbers that can be used to characterize the 
problem of PDR structure. Such estimates, and the ac- 
companying physical insights, would allow modeling of 
clouds over a wide range of galactic environments with- 
out the need for a complex and numerically costly radia- 
tive transfer calculation to cover each case. 

Our goal in this work is to revisit the problem of de- 
termining the size of a PDR in a finite cloud embedded 
in a multi-dimensional radiation field, and to derive an- 
alytic approximations for the structure of a PDR that 
will yield gross yet observable quantities such as the to- 
tal atomic hydrogen column around a molecular region 
and the fraction of a cloud's volume in the atomic and 
molecular phases. As part of this work we determine the 
important dimensionless numbers that characterize the 
problem, and we provide a rough classification of PDRs 
based on them. The model we develop is capable of span- 
ning the range from galaxies where the gas in atomic- 
molecular complexes is predominantly molecular and a 
slab treatment is appropriate to dwarf galaxies where 
only a tiny fraction of the ISM resides in the molecular 
phase. In future work (Krumholz, McKee, & Tumlinson, 
in preparation), we provide a more detailed application 
of the results derived here to the problem of determining 
the atomic to molecular ratios in galaxies. Before moving 
on, we do note that our focus on an analytic solution with 
a multi-dimensional radiation field, characterized by a 
few dimensionless numbers, has a price: our approach to 
the chemical and thermal physics of PDRs is significantly 
simpler than much previous work. We do not account for 
factors such as the temperature-dependence of rate coef- 
ficients or H2 dissociation by cosmic rays. Our work is 
therefore less suited to making detailed predictions of the 
structures of individual PDRs than it is to making pre- 
dictions for galactic-scale trends in atomic and molecular 
content. 

We approach the problem of finite clouds by idealizing 
to the case of a spherical cloud embedded in an isotropic 
radiation field, since this allows us to explore the effects 
of finite cloud size and curvature while at the same time 
keeping the problem simple enough to admit an approxi- 
mate analytic solution. Our approach is as follows: in §[5] 
we state the formal problem and introduce some physi- 
cal approximations that are independent of geometry. In 



§ [3] we derive an approximate analytic solution to the 
one-dimensional semi-infinite slab case, which allows us 
to demonstrate the underlying physical principles of our 
approach. In this section we also compare to a grid of nu- 
merical solutions and show that our approach produces 
good agreement. Then in § [4] we extend our approach 
to handle the case of a spherical cloud embedded in an 
isotropic radiation field. Finally we summarize and draw 
conclusions in § [H 

2. THE FORMAL PROBLEM 

Consider a region of hydrogen gas where the number 
density of hydrogen nuclei is n, mixed with dust which 
has a cross section a^^ per H nucleus to radiation of fre- 
quency v. The hydrogen is a mix of atoms and molecules, 
with a fraction /hi of the nuclei in the form of Hi and a 
fraction /h 2 = 1 — /hi in the form of H2. We consider 
frequencies v that fall within the Lyman- Werner (LW) 
band from v\ = c/1120 A to V2 = c/912 A, such that 
photons of that frequency can be resonantly absorbed by 
hydrogen molecules. 

The equation of radiative transfer for a beam of radi- 
ation in direction e passing through this gas is 



e • VI V = —n I 7t/h 2 cth 2 ,!/ + o&.y ) l v 



(1) 



where l v is the radiation intensity at frequency v and 
CTH 2 ,iy is the cross section for absorption of radiation at 
frequency v by a molecule of hydrogen. The value of 
(Th 2 .i; may change with position as the fraction of H2 
molecules in different quantum states changes. The total 
fraction of the gas in the molecular phase is determined 
by the balance between the rate of H2 formation and 
dissociation, 



fmn 2 K = ^-n J dfl J dv ^cr H2 ,^/ 



(2) 



where 1Z is the rate coefficient for formation of H 2 
molecules on dust grain surfaces and /diss,!/ is the frac- 
tion of absorptions at frequency v that yield dissociation 
of the H2 molecule rather than decay back to a bound 
state. 

Note that we do not include a source term in the trans- 
fer equation |T]), because although most FUV photons 
absorbed by H2 molecules do decay through a vibrational 
ladder via photon emission, the photons released in this 
process do not fall into the LW band. Thus, the trans- 
fer equation we have written is only valid for frequencies 
in the LW band. We have also neglected scattering of 
FUV photons by dust grains. Since scatt ering is highly 
forwa rd-peaked at FUV wavelengths (e.g. iRoberge et alj 
119811) . this approximation is reasonable as long as we 
take Ud.u to be the absorption cross section, not the total 
cross section. We have also omitted other H 2 dissociation 
mechanisms other than LW photons, such as cosmic ray 
collisions and chemical reactions. These are significant 
only in nearly fully molecular regions where there are no 
significant numbers of LW photons present. 

Equations fl} and ((2J, together with the atomic and 
dust physics that specify cth 2 ,i> and a^ v and a boundary 
condition that specifies I v on all rays entering the surface 
of a cloud, fully determine /„ and /h 2 at all positions. 
We cannot solve them exactly, but we can obtain an ap- 
proximation that exposes the basic physical outlines of 
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the solution. We b egin by making two stand ard approxi- 
mations, following iDraine &; Bertoldl (|1996l ). to simplify 
the atomic physics. 

First, it is convenient to simplify the transfer equation 
(fTj) by dividing by hv to transform from intensity to pho- 
ton number, and then by integrating over frequency in 
the LW band. In so doing we can exploit the fact that 
for realistic dust <Jd.v is nearly independe nt of frequency 
in the LW band (jDraine fc B crtoldi 1996) to replace er^ 
with a constant value <7d- Doing so gives 



e-Vf 



-ntJdl* 



1 , 
-™/h 2 



duan 2 , v It, (3) 



where /* = I v /(hv) is the photon number intensity, i.e. 
the number of photons per unit time per unit area per 
unit solid angle per unit frequency that cross a given sur- 
face, and I* = dv I* is the photon number intensity 
integrated over the LW band. 

Second, we note that /diss,^ varies only weakly when in- 
tegrated over frequ e ncy a nd over position within a PDR. 
IDraine fe~ Bcrtoldi (1996|) show that over the width of 
a PDR it stays roughly within the range 0.1 — 0.2. Its 
value in f ree space depends on th e assumed radiation 
spectrum; IDraine &: BertoldH (|1996h fin d /diss = 0-12 in 
free sp ace for their fiducial choice, while iBrowning et alJ 
(2003) suggest /diss = 0.11 as a typical value. For sim- 
plicity we adopt a constant value fdiss,u — /diss = 0.1 
and take this constant out of the integral, reducing the 
dissociation equation to 

f H in 2 U = f J^l n J dV. £ diyi;a S2 , u . (4) 

It is convenient at this point to produce a combined 
transfer-dissocation equation from Q and (j4|) . If we in- 
tegrate equation {3| over solid angle dfl, we obtain 

V • F* = -na d cE* - ^f B2 n J dfl J J dvau 2 , v I* v , (5) 

(6) 



where 



E* = - I ,1V. I 



F* 



dV. el* 



(7) 



are the photon number density and photon number flux 
integrated over the LW band, respectively. We can then 
use equation (01 to substitute for the last term, yielding 
the combined transfer-dissociation equation 



V • F* = -na d cE* - 



fmn 2 7l 



diss 



(8) 



3. SOLUTION IN ONE DIMENSION 



3.1. Analytic Solution 

We start by giving an approximate analytic solution 
to this formal problem for unidirectional radiation im- 
pinging on a one-dimensional semi-infinite slab in order 
to illustrate the physical principles behind our approach. 
Consider a region of gas of density n filling the half-space 
z > 0, subjected to a dissociating radiation field of pho- 
ton number intensity /* = 47rjQ<5(|e — z\) that fills the 
half-space z < 0, where Jq is the angle-averaged inten- 
sity in free space. The corresponding free-space photon 



number density is Eq = An Jq / c, and the magnitude of 
the free-space photon flux is Fq = cEq. For simplicity 
we neglect the (relatively weak) temperature-dependence 
of ft. 

Since the radiation intensity everywhere at all z re- 
mains proportional to <5(|e — z\), it immediately follows 
that 

F* = F*z = cE*z, (9) 

at all points, and the combined transfer-dissociation 
equation reduces to 

dF* „ fmn 2 Tl 

— — = -na d F - 
dz 



(10) 



/diss 

subject to the boundary condition that F* = Fq at 
z = 0. Since numerical calculations show that the tran- 
sition from predominantly atomic gas to predominantly 
molecular gas in a PDR generally occurs in thin band 
bounded by much larger regions where the gas is either 
predominantly atomic or predominantly molecular, we 
can obtain a good approximation to the exact solution 
by treating /hi as having a constant value near unity over 
the bulk of the PDR, and then dropping to zero as a step 
function once the fully molecular surface is reached. For 
constant /hi = 1, we can non-dimensionalize equation 

CEo]) to 

dr x 

where T — F* / (Fq) is the fraction of the free-space flux 
remaining, r = nodZ is the dust optical depth from the 
slab surface, and 



X 



/dissCrdC-Ep 

nR 



Equation (|11[) has the exact solution 



X 



-(t-thi) 



(12) 



(13) 



where 



thi = ln(l +x) 



(14) 

is the depth at which the flux goes to zero, which we take 
to be the optical depth through the Hi region. Of course 
in reality the flux should never go to zero exactly. That 
is does in our solution is an artifact of our choice to treat 
/hi as constant. Nonetheless, since the transition from 
/hi w 1 to /hi w is sharp, z Ha = r H i/(ncr d ) should 
be a good approximation of the depth at which the gas 
becomes predominantly molecular. 

The dimcnsionless parameter x/ fm is the ratio of the 
two terms on the right hand side of equation (fT0|) with F* 
set equal to its value AitIq at the slab edge. This makes 
its physical meaning clear: x/fwi represents the ratio 
of the absorption rate of LW photons by dust grains to 
the absorption rate by H 2 molecules for a parcel of gas 
exposed to the unattenuated free space radiation field. 
If the gas at the edge of free space is predominantly 
atomic, as is the case for example at the edge of an 
atomic-molecular complex, then /hi ~ 1 and this ratio is 
simply given by x- F° r X > 1 absorptions by dust grains 
dominate, while for x < 1 absorptions by H2 molecules 
dominate. 

For a giant molecular cloud in the Milky Way and its 
outer atomic envelope, typical values of the number den- 
sity, dust cross section, and H2 formation rate coefficient 
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are n ~ 30 cm" 3 , a d ~ lO"^ 1 cm z and Tl ~ 3xl0~ u cm" 
s _1 (|Draine fc BertoldilH99l . Using the [Draini jT97l) 
functional form for the local FUV radiation energy den- 
sity as a function of wavelength, 

14 A3 5 (31.016A|j 

-3 



A£ A = 6.84 x 10 

erg cm 



49.913A 3 



19.897) 
(15) 

the free-space photon number density from 912 — 1120 
A is £q = 7.5 x 10~ 4 cm~ 3 . For an H 2 molecule in the 
ground state, this corresponds to a free-space dissocia- 
tion rate of 3.24 x 10~ n s _1 . (In principle for a slab 
computation we should divide the observed value of Eq 
by two to account for the fact that one can only see half 
the sky at the surface of an opaque cloud, but we do 
not do so here because in § |4] we will account for this 
effect self-consistently.) Thus, for Milky Way conditions 
not near a local strong source of FUV, x °f order a few 
might be typical. Thus, in the Milky Way dust shielding 
is marginally significant in determining the structure of 
atomic-molecular complexes. 

3.2. The Two-Zone Approximation 

We can integrate the transfer-dissociation equation (JSJ) 
directly in one dimension because, due to the constant 
angular distribution of the radiation, there is a trivial re- 
lationship between E* and F*. In multiple dimensions, 
however, there is no simple relationship between the two, 
because the angular distribution of the radiation inten- 
sity is not constant with position inside a cloud. To over- 
come this problem, we adopt what we call the two-zone 
approximation. When the photon number density E* is 
large, the first term on the right hand side of equation 
((U, representing absorptions of photons by dust grains, is 
much larger than the second term, representing absorp- 
tions by H2 molecules. This makes intuitive sense: in 
regions where there many dissociating photons present, 
the molecular fraction will be very low, so there will be 
few H2 molecules available to absorb LW photons and 
most photons will be absorbed by dust. In regions where 
E* is small, the number of molecules will increase, and 
for any appreciable number of molecules these will dom- 
inate the absorption rate. 

In the two-zone approximation, we divide the cloud 
into a zone where dust absorption dominates and a zone 
where molecular absorption dominates. In the dust- 
dominated region we drop the molecular absorption term 
in the radiative transfer or transfer-dissociation equa- 
tions (equations [3] or [8|) , and approximate the opacity as 
having a constant value nod- In the zone where molecular 
absorptions dominate, we drop the dust absorption term 
in equation j8J) and approximate the molecular absorp- 
tion term by (ftn 2 lZ/ /diss, where (ft > 1 is a constant of 
order unity, whose precise value we determine below, that 
we include to account for the fact that some LW photons 
will be absorbed by dust grains even in the molecular- 
dominated region. We define a boundary between these 
two zones by the condition that the dust and molecular 
absorption terms be equal, which is satisfied when 



E* 
E* 



X 



1 



(16) 



where we have set /hi = 1 because at the point of 
equality the molecular fraction is <C 1, and for conve- 



nience we have defined the modified dust to molecu- 
lar absorption ratio ift — x/4>- With this approxima- 
tion, the one-dimensional non-dimcnsionalized transfcr- 
dissocation equation becomes 



dT 
~dr 



F, F>l/ift 
l/ift, F < 1/ift 



(17) 



We shall see in §H]how the two-zone approximation en- 
ables us to solve the problem in the spherical case. First, 
though, we examine the solution in the one-dimensional 
case. If ip < 1, then F < 1/ift is satisfied everywhere and 
equation (flTf has the trivial solution 



F : 



4> 



(18) 



The flux goes to zero at a depth thi = ift. If ift > 1, the 
solution is 



F 



with 



e , t < r d 

V'~ 1 (' r HI - t)/(t H I - T d ), T > T d 

t ( \ = In ift 
thi = 1 + hi ift. 



(19) 



(20) 
(21) 



Here Td represents the dust depth into the slab at which 
the absorption begins to be dominated by H2 molecules, 
while thi is the optical depth where we expect a transi- 
tion from mostly atomic to mostly molecular gas. Com- 
bining the two cases, we have 



thi = 



ift, ift < 1 

1 + In ift, ift > 1 



(22) 



We now turn to the question of determining the con- 
stant (ft. Physically, we expect to have <j> — * 1 for x *C 1; 
because in that case dust absorptions contribute negligi- 
bly throughout the cloud. For x> 1 we expect to have 
(ft to asymptote to a value greater than unity, account- 
ing for the contribution of dust to absorptions even in 
the molecular-dominated region. A comparison of the 
limiting behavior of the analytic solution (jT4j) with the 
two-zone approximation (|22|l confirms this physical ar- 
gument, and suggests that the appropriate limiting be- 
havior is (ft — * 1 as x —> and (ft — ► e as x ~ * 00 ■ We 
therefore adopt 

- 2 ' 5 + * e (23) 



2.5 



■ X 



which has the correct limiting behavior, and where the 
value 2.5 is chosen to optimize agreement between the 
two-zone approximation and the analytic solution in the 
intermediate x region. 

3.3. Comparison to Numerical Calculations 

Before using the two-zone approximation to compute 
the case of a finite cloud with an isotropic radiation field, 
we check its accuracy for the one-dimensional case by 
comparing with detail ed numerical calculations using the 
iBrowning et alj (|2003t ) H2 formation and radiative trans- 
fer code. We refer readers to that paper for a full de- 
scription of the physics included in this calculation, but a 
brief summary is that the code numerically integrates the 
frequency-dependent equation of radiative transfer for a 
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unidirectional beam of radiation incident on an isother- 
mal, constant-density slab of gas mixed with dust. The 
transfer equation is coupled to a statistical equilibrium 
calculation that determines the populations of Hi atoms 
and a large number of rotational and vibrational levels 
of the H2 molecule that are excited by LW band photons 
in each computational cell. The output of this calcula- 
tion is the fraction of H nuclei in molecules as a function 
of depth within the cloud. The code we use her e dif- 
fers from that described in 1 Browning et al.l (j2003j) only 
in that the earlier version accounted for absorptions of 
LW photons by dust grains by mo difying the photodisso- 
ciation rate using the method of Ivan Dishoeck &: Blackl 
(1986]), whereas the version we use here computes radia- 
tion attenuation by dust grains directly from the radia- 
tive transfer equation. 

For the models we present here we use a density and 
temperature of n = 5000 cm -3 and T = 90 K. These 
values are chosen purely for computational convenience, 
and have no significant impact on the results. The in- 
cident radiation field is a unidirectional beam of pho- 
tons uniformly distributed in frequency over the wave- 
length range 912 — 1120 A. The frequency-dependent 
photon flux in this beam is F*, so Eq — F*(v% — v\)/c 
and Jq = F*{v2 — vi)/(An). We adopt a dust extinc- 
tion c urve following the functional form of Cardclli et al. 
(1989), scaled to give a dust cross section per H nucleus 
at 1000 A of a = &d,MwZ', where Z' is the metallicity 
relative to solar and we take (7d,MW = 6.0 x 10~ 22 cm 2 or 
2.0 x 10~ 21 cm 2 to be two fiducial dust opacities for the 
Milky Way. These two values of <7d,MW correspond to the 
est imated attenuation cross sections at 1000 A estimated 
bv lDraine fc Bertoldil ([1996) for dense and diffuse clouds 
Milky Way, respectively. We adopt a rate coefficient for 
H2 formation on grain surfaces of 1Z = Kmw^' with 
K M w = 3 x 1Q- 17 Z' cm" 3 s _1 as our fiducial Milky Way 
value ( Wol fire et al.ll2008h . We do our computations for a 
grid of models running from F* = 10~ 7 — 10~ 3 photons 
cm~ 2 s _1 Hz -1 in steps of nine steps of 0.5 dex, from 
Z' = 10~ 2 — lO 5 in six steps of 0.5 dex, and for the two 
values of (Td,MW mentioned above. These values F* are 
significantly higher than are typical in the Milky Way, 
but are chosen so that, in conjunction with our choice of 
n, the ratio E*/n that appears in x is within the typical 
Milky Way range. With this parameterization 

X = 0.75(CTd,MW,-2l/7?-Mw)(^0 l n 2) (24) 
= 4.07CT d , MW ,-2lF*-5 (25) 

where Od,MW,-2i = tfd,Mw/10 -21 cm 2 is the Milky Way 
1000 A dust absorption opacity normalized to 10 -21 cm 2 , 
Eq = Eq/7.5 x 10~ 4 cm~ 3 is the free-space dissociat- 
ing photon number density normalized to the Milky Way 
value, ?i2 = n/100 cm~ 3 is the number density of hydro- 
gen nuclei in units of 100 cm~ 3 , and F*_ 5 = F*/10~ 5 

photons cm~ 2 s _1 Hz -1 . Thus the calculation covers a 
broad range of parameters from strongly dust-dominated 
to strongly molecular-dominated, thereby bracketing the 
real Milky Way value of x ~ 1- Note that x is inde- 
pendent of Z' because for the parameterization we have 
chosen the Z'-dependences of Cd and 1Z cancel. Since we 
predict that the dust optical depth through the PDR, 
Thi = nadZn 2 , depends only on x, and x in turn depends 
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Fig. 1. — The plots show /h 2 versus dust optical depth r = na^z 
for our numerical radiative transfer calculations with log Z' = —1.5 
{dashed lines), logZ' = —0.5 (solid lines), and logZ' = 0.5 (dot- 
dashed lines). The gray vertical lines indicate the optical depth of 
the transition to fully m olec ular as calculated with the two-zone 
approximation, equation il2_'l) . Each cluster of three curves plus a 
vertical line indicating a prediction corresponds to a radiation flux 
logF* = —7, —5, or —3, as indicated. The two panels are for the 
cases <r d mw ,—21 = 0-6 and 2.0, as indicated. 

only on the ratio <Jd/7l, thi should be independent of 
Z. Since we use a range of 10 2 5 in Z' , our numerical 
calculations represent a strong test of this prediction. 

Figure [T] shows /h 2 versus depth within a cloud as com- 
puted numerically for a sample of our input parameters, 
overlaid with the corresponding locations of the atomic 
to molecular transition as calculated via the two-zone 
approximation. As the figure shows, the two-zone ap- 
proximation does a very good job reproducing the loca- 
tion of this transition over an extremely broad range of 
parameters. 

To quantify the quality of the approximation, we must 
define a fiducial measure for the depth of the H2 region 
in the numerical calculations, since in these runs /hi ap- 
proaches but never reaches unity. The most reasonable 
measure is 

/■OO 

N m = / dzfnm, (26) 
Jo 

the total Hi column integrated through the cloud. Since 
the radiation field is attenuated exponentially or faster, 
and /hi is proportional to radiation intensity in the re- 
gion where /h 2 ~ 1, this integral is guaranteed to con- 
verge. In the limit where the transition from HI to H2 
is sharp it approaches the total gas column up to the 
transition point. In practice we cannot continue the nu- 
merical integration to z = 00, so we truncate the integral 
at the value of z where /hi = 5 x 10~ 3 ; using /hi = 10~ 2 
instead changes the value by less than 8% for all our 
runs, and by less than 2% for all runs with \ > 0.1, so 
our evaluation of the integral should be accurate to this 
level. 

We plot the dust opacity through this hydrogen col- 
umn, iVniCTd, and the corresponding value thi predicted 
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Fig. 2. — The upper panel shows the dust opacity Affjicd through 
the numerically-determined Hi column (various symbols) as a func- 
tion of x- The values of Z' and ct^mw — 21 f° r each calculation are 
indicated by the plot symbol. This is compared to the optical depth 
tjji computed from the two-zone approximation (equation [22] thin 
solid line) and computed using the analytic solution (equation 1141 
thick dashed line). The lower panel shows the error in the two-zone 
approximation, defined as Error = iVm^d/THl — lj the dotted line 
indicates zero error. 

by the two-zone approximation, in Figure [2j As the 
figure shows, the two-zone approximation recovers the 
numerically-computed Hi column to better than 50% ac- 
curacy over almost a five-decade range in \- The error in 
the two-zone approximation is generally comparable to 
or smaller than the spread between models with different 
dust opacities but the same value of x- 

The error in our approximation is largest at small x, 
and examination of Figure [1] suggests the reason why: by 
evaluating the equations with /hi = 1 inside the PDR, we 
have assumed that the transition from atomic to molecu- 
lar is sharp. This is true for \ ~ 1 or greater, but begins 
to fail for \ <C 1. In our runs with x > 1j typically 95% 
of the gas is atomic in the region where /hi > 0.5; and 
even at the depth where /hi drops to 5 x 10 -3 more than 
half the gas column above that point is atomic. This in- 
dicates a very sharp atomic-molecular transition, so our 
approximation that /hi = 1 until the gas is almost en- 
tirely molecular is a good one. For \ ~ 0.01, on the other 
hand, roughly 80% of the gas at /hi < 0.5 is atomic, and 
Hi contributes only 10% of the total gas column above 
/hi = 5 x 10~ 3 . The transition from atomic to molecular 
is therefore much more gradual, and our accuracy suffers 
as a result. 

Nonetheless, we note that x « 1 does not appear to 
be physically realized in normal galactic environments. 
For Milky Way molecular clouds \ ~ 1 or greater, and 
reducing x to 0.01 would require some combination of 
reducing the ISRF and increasing the atomic gas den- 
sity by a factor of 100. Such a combination of very high 
atomic ISM density and very low radiation field is gener- 
ally not observed. We conclude that, for realistic physical 




Fig. 3. — Illustration of the two-zone approximation in spherical 
geometry. 

parameters, and given that these parameters (such as Co 
and 1Z) a re themselves unce rtain at the factor of a few 
level (e.g. IWolfire et al.ll2008l ). the error in the two-zone 
approximation is unlikely to be the dominant one. 

4. SOLUTION FOR SPHERICAL CLOUDS 

We now extend the two-zone approximation to a spher- 
ical cloud of radius R embedded in a uniform, isotropic 
radiation field of angle-averaged intensity Jq . (Note that 
this radiation field has the same LW photon number den- 
sity as the unidirectional radiation field considered in §[3l 
so it gives the same dissociation rate in free space.) Fig- 
ure [3] illustrates the basic geometry of the problem and 
our approximation: we consider the dust-dominated re- 
gion to run from radius r = to r = R, and the molec- 
ular self-shielding region to run from r = rn 2 to r — r<j. 
For convenience we introduce the dimensionless position 
variables x — r/R and y = 1 — x, and we define the dust 
optical depths from the surface to and to th 2 as Td 
and 7~hi, respectively. 

In § 14.11 and § 14.21 we develop the basic equations that 
describe the two-zone approximation for clouds with and 
without dust opacity-dominated envelopes. We then ex- 
plore three limiting cases of these equations. We consider 
the behavior at the boundary between the presence and 
absence of a dust-dominated zone in § 14. 3[ and we ex- 
plore several interesting limits in § 14.41 We then give a 
numerical solution and an analytic approximation to it 
in § 14.51 In § 14.61 we compare our solution for a finite 
cloud to the standard slab approximation, to determine 
when the slab approximation is valid and when it fails. 
In § I4.7[ we address the level of uncertainty introduced 
by the approximations we make in the spherical case. 
Finally, in § 14.81 we present some example calculations 
using our analytic approximation. 

4.1. Clouds with Dust- Dominated Zones 

First consider the case where ip is large enough so 
that there is a dust-dominated zone in the outer part 
of the cloud where molecular self-shielding is negligible, 
i.e. rd < R. The transfer equation in this region becomes 

e- VI* = -na d r, (27) 

which for rays originating at the cloud surface and stay- 
ing entirely within the dust-dominated region has the 
trivial solution 



I*{x,n) = exp(-r R ^)J *, 



(28) 
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where tr = na^R is the center-to-edge dust optical depth 
of the cloud, 



£ = — X 2 + X 2 /l 2 — Xflj , 



(29) 



is the distance, normalized to the cloud radius, from ra- 
dius r to the cloud surface on a ray that makes an an- 
gle 9 relative to the radial vector (see Figure [3]), and 
\i = cos 9 = — e ■ f. This solution applies for /i > 0. On 
the other hand, if fi < /ih 2 , 



/^H 2 



xn 2 

Xd 



(30) 



then the ray passes through a part of the cloud that is 
fully molecular. Since the fully molecular region will be 
extremely opaque, to good approximation along these 
rays I*(x,^) = 0. Finally, rays for which > // > /in 2 
pass through the region where the opacity is dominated 
by molecules rather than dust, but where the gas is not 
yet fully molecular. Again, we use the approximation 
that the transition from a low molecular fraction to fully 
molecular is sharp, so that over most of this region the 
molecular fraction is not vastly larger than it is at the 
region's edge. This is consistent with numerical solu- 
tions of the problem (e.g. Figure [I}, which show that 
the atomic fraction /hi is nearly constant through the 
bulk of the PDR, and rises from its free space value to 
unity over a small region. This means that the opacity 
is not much greater than its value of nod at the outer 
edge of the molecule-dominated region. We therefore 
approximate that the optical depth along these rays is 
the same as for those in the dust-dominated region, tr,£. 
This approximation is perhaps the least certain part of 
our calculation, and we quantify the level of uncertainty 
that it produces in § 14.71 Combining these three regions, 
we have an approximate intensity 



exp(-T R £)J *, > /^h 2 

0, [i < /1H 2 



(31) 



E*(x d ) 
EX 



The location of Xd, the crossover point from dust- to 
molecule-dominated absorption, is defined by the condi- 
tion that E*/Eq = 1/ip (equation [TH)) . For convenience 
we define the function 770 by 

1 f 1 

??o(a;d,^H 2 ;TR) = - / dfj, exp (-r R £ d ) , 

(32) 

where £d = £(^d)- We show that 770 can be evaluated in 
terms of exponential integrals in Appendix [AJ The first 
equation for the two-zone approximation is therefore 

Vt)(xd,xn 2 ;T R ) = (33) 

Inside Xd, we drop the dust opacity term, so that the 
combined transfer-dissociation equation ([8]) becomes 

where we have again assumed that /hi ~ 1 outside the 
fully molecular region, and for convenience we have in- 
verted the sign by defining F* = —F*r. The solution 
is 

'Xu 2 



= 3~V> 



1 - 



where we have chosen the constant of integration by re- 
quiring that F* — at x — xn 2 . To determine a;n 2 from 
the boundary conditions, however, we must determine 
the flux at some other location. Thus we evaluate the 
flux F* at Xd by integrating the intensity over solid an- 
gle using equation (|3"Tj) . For convenience we define the 
function i]i by 

F*(x d ) 



F* 



1 f 1 

r)x(xd,xn 2 ]TR) = - / dfi fi exp ( Tf^d) . 

1 J Uw„ 



■ m (36) 
As with 770, in Appendix IAI we evaluate 771 in terms of 

exponential integrals. Combining (|36p with (|35| gives an 

implicit equation for xn 2 : 



T]i(xd,xii 2 ;T R ) 



TRXd 

3-0 



1 



2-H2 
Xd 



(37) 



> 



(35) 



Together, equations (|3"3"|) and (|3T|) constitute two equa- 
tions in the two unknowns Xd and xh 2 , and thus fully 
determine the location of the transition from predomi- 
nantly atomic to predominantly molecular in the two- 
zone approximation. 

4.2. Clouds without Dust- Dominated Zones 

Now consider the case where t/> is small enough so that 
there is only one zone, because even gas at the edge 
of the cloud is sufficiently molecular for absorptions by 
molecules to outnumber those by dust grains. In this case 
equation (|35[) applies throughout the cloud, so we must 
fix xh 2 directly from the boundary conditions. To do 
so we need to know the flux F*(l) at the cloud surface. 
This is not simply Fq = cEq as in the case of a uni- 
directional radiation field; in free space for an isotropic 
radiation field F* vanishes, and _F*(1) is non-zero only 
because rays passing through the cloud do not carry the 
same intensity as rays that do not pass through it, pre- 
venting the integral over angle from vanishing. Thus for a 
sufficiently transparent cloud, F*(l) approaches zero, its 
value in free space. Conversely, at the surface of a cloud 
that is opaque and extremely large, F*(l) = Fq/4. The 
factor of 1/4 relative to the unidirectional case arises 
because because half the solid angle is blocked by an 
opaque object (providing one factor of 1/2), and because 
in the part of the sky that is not blocked the radiation is 
isotropic, and one must average over all the directions in 
which photons are traveling to find the fraction of that 
motion in the —f direction (providing another factor of 
1/2). 

The problem of determining the intensity is exactly 
the same as in § 14.11 At the surface of the cloud rays at 
angles /1 > do not pass through the cloud and there- 
fore contribute the unattenuated free space intensity Iq . 
Those with [i < /ih 2 = — (1 — ;e h 2 ) 1 P ass through the 
fully molecular region and therefore contribute zero in- 
tensity. For rays at angles /in 2 < /i < 0, we make the 
same approximation as in § 14. li that the molecular ab- 
sorption rate per unit distance that a photon travels is 
roughly constant until one approaches the sharp tran- 
sition from atomic to molecular. Thus the molecular 
opacity-dominated part of the PDR has a constant ef- 
fective opacity, which we can determine by computing 
its value at the cloud surface. For convenience we char- 
acterize this opacity via an effective cross-section per H 



nucleus o~ e . By examining the transfer-dissociation equa- 
tion ((HJ, it is clear that for /hi ~ 1 this opacity is 

drnlZ 



f diss cE*(iy 



(38) 



With this approximation the transfer equation through 
the region outside where the gas becomes fully molecu- 
lar is simply equation (|27|) with <7d replaced by cr e , and 
the solution for the intensity along each ray is given by 
equation (|31j) with r R replaced by r cR = {o~ c / o~&)tr. The 
photon number density and flux at the cloud surface are 
therefore given by 



E*(l) 
F*(l) 



= r/o(l,x H2 ;r cR ) 
= 771(1, xh 2 ;t cR ). 



(39) 
(40) 



With these arguments, evaluating equations (|A11|) and 
(IA12P shows that 770 and 771 reduce to 



r7 (l, ih 2 ;t c r) = 
771(1, XH 2 ;r e R) = 



4r^ + 2 
(1 - 2MH 2 T eR ,)e 2 ^ r « - 1 



(41) 



r-(42) 



Using these values of E*(l) and F*(l) in equations 
and (l35l), we find 



1 



2/.iH 2 r cR 



4t cR , 

(1 - 2 MH2 r eR )e 2 ^^oH 
8r e 2 R 



3ip 



1 



r 3 
f H 2 



(43) 
.(44) 



We therefore again have two equations in two unknowns, 
with the unknowns in this case being r oR and xr 2 . Note 
that in this case /ljh 2 can be isolated in the first equation, 
to give 



1 



Mh 2 



2t, 



cR 



In 1 + 2r oR 



4tr\ 

i/> J 



Together with the relation between /jh 2 and ih 2 (equa- 
tion I30p , this reduces the problem to a single non-linear 
equation, which is convenient for numerical solution. 
However, the two-equation form is more convenient for 
an analytic approach. 

At this point it is worth making a few remarks about 
the behavior of equations (|4"3"| and (|4"4")l . First, equation 
(|43|) implies that t r /(t/jt cR ) > 1/2, so the argument of 
the logarithm in equation (|45|) is always less than unity 
and /Ltn 2 is negative. Second, in all of these equations 
tr and ip appear only through the combination tr/iP, 
so values of Xn 2 and T eR must be constant on lines of 
constant tr/tJj. Finally, note that in § 14.51 we give an 
approximate analytic solution to equations (|43[) and (|4"4"|) . 

4.3. The Dust- Dominated Zone Boundary 

We first identify the boundary between the presence 
and absence of a dust-dominated region. In the case of 
a perfectly beamed radiation field impinging on a semi- 
infinite planar slab, which we treated in §[31 this is ij) = 1. 
The result is more complex in the case of an isotropic 
radiation field and a cloud of finite size. The bound- 
ary between the two cases is defined by the condition 



Dust zone 
No dust zone 



Fig. 4. — The thin black curves show the boundaries in the 
( r Ri VO-plane at which x d = 1, t c r = tr, and at which xjj 2 = 0, as 
indicated by the text accompanying each curve. The thick dashed 
curve shows our approximation to the x<j = 1 curve, equation H47I I. 
The asterisks along the curve for x& = 1 mark the points at which 
xh, = 0.1, 0.3, 0.5, 0.7, and 0.9, as indicated. 

that Xd = 1 or r cR = r R , i.e. that the cross-over be- 
tween dust-dominated and molecular-dominated absorp- 
tions occur at the cloud surface, or equivalently that the 
dust and molecular effective opacities at the cloud sur- 
face are equal. It is immediately obvious that equations 
and (|3"T|) become identical to equations (|4"3"]l and 
in this limit. If we set r eR = r R , then equations (|45|) 
and ([44]) define a curve in the (r R , ?/')-plane that cor- 
responds to the point where dust-dominated layer dis- 
appears. Above the curve the radiation field is strong 
enough so that the outer part of the cloud is dust opacity- 
dominated, while below it molecular opacity dominates 
throughout. Along the bounding curve, 



• T H 2 = \ 1 - 



1 



In 1 



-1 2 



(46) 



(45) ® 



We show the curve at which the dust-dominated layer 
vanishes, and the value of xr 2 along the curve, in Figure 



That along this curve ip — > 1 for r R <g; 1 and tp — > 2 for 
tr 1 makes intuitive sense. If tr 1, the dust optical 
depth through the cloud is tiny and so the radiation field 
has its unattenuated, free-space value regardless of posi- 
tion in the cloud. Thus the condition that dust shielding 
and molecular shielding contribute equally at the cloud 
surface (xd = 1) can only be fulfilled is they are equal 
or nearly so in free space, which is simply a requirement 
that ip = 1. Similarly, if t r 3> 1 then the cloud is effec- 
tively a semi-infinite slab whose curvature is negligible. 
In this case the radiation field at the cloud surface only 
contains contributions from rays with /1 > 0, i.e. those 
that never pass through the cloud; rays with fi < are 
infinitely attenuated. Thus, the radiation field has ex- 
actly half its free-space value, the molecular fraction has 
double its free-space value, and the condition that dust 
shielding and molecular shielding are equal reduces to 
the requirement that ijj = 2. 

These considerations suggest that a function that in- 
terpolates between these two limiting behaviors is likely 
to produce a good approximation. Numerical experimen- 
tation shows that the curve 



1.4 + 2tr 

1.4 + Tr 



(47) 



9 



reproduces the true value of ip along the curve Xd = 
1, T G R = tr, to better than 3% for all tr. We show 
this approximate solution with the thick dashed curve in 
Figure [U 

4.4. Limiting Cases 

We can better understand the behavior of PDRs in 
finite clouds by exploring several limiting cases of our 
equations, corresponding to clouds that are very large or 
very small, and to radiation fields that are very strong 
or very weak. 

Case 1: Strong Radiation Fields. The first limit we 
consider is one in which the radiation field is so strong 
that there is no fully molecular core, so xr 2 = 0. It is 
easy to verify that if there is no dust-dominated region, 
so equations (|43|) and (|44p apply, then there are no fi- 
nite values of ip and tr such that xn 2 = 0. (However, 
see § 14.71 where we show that this behavior is probably 
not physical.) On the other hand, if there is a dust- 
dominated region, equations (|33|) and (|37|) apply and 
x n 2 = can be reached at finite ip and tr. This be- 
comes clear if we note that xn 2 = implies //h 2 = — 1 
(following equation l30|) . and equation (f37|) then admits 
the solution Xd = 7/i(xd, 0; tr) = 0. Since Xd = 0, it 
immediately follows that £d = 1 and 770(0, 0;tr) = e _TR . 
Equation ([33]) then gives ip = e TR . The physical meaning 
of this solution is that ip — e TR is the critical curve along 
which xh 2 = Xd = 0; at this value of ip or larger the 
radiation field is too intense for a fully molecular core 
to exist. We plot the critical curve in Figure HI In Ap- 
pendix[B]we solve equations ([33]) and (|37|) perturbatively 
in the vicinity of the critical curve, and show that near 
the strong radiation boundary the solution is 



Xd : 



xh 2 



G 



tr(t r + 2) \ iIj 
1536 



1 



25tr(t r + 2f 



1/4 



1/2 



\ 5, 



(48) 
(49) 



This solution obviously only applies for e TR > ip. 

Case 2: Small Clouds. Before analyzing this case, we 
warn that in ? 14. 71 we show that our solution in this case 
should be regarded as giving an upper limit on the molec- 
ular fraction rather than a direct estimate. However it 
is still useful to consider this case, both in order to de- 
rive upper limits and to provide expressions that can be 
incorporated into approximations in parts of parameter 
space where our method does provide estimates rather 
than upper limits. We have shown that for ip > 1, there 
is a finite value of tr at which the fully molecular core 
vanishes, and conversely that if ip < 1 there is no finite 
tr for which xn 2 — 0. However, one can easily verify 
that when there is no dust-dominated zone (ip < 1) and 
equations (|43| and (|44|) apply, ih 2 — ► as tr — > 0. Phys- 
ically, this corresponds to the case of a cloud that is so 
small that its dust is optically thin to LW photons. In 
Appendix [Cl we show that in this limit the solution may 
be approximated by 



t c r ■ 



Xn 2 



TR 

IP ' 



2ip 2 4ip 3 

2VM 



1± 

60ip 4 



+ 



VZip 5^ 



(50) 
(51) 



Note that by definition t cR /tr = [Eq / E* (1)] / 'ip , so 
equation (fSTJ)) is effectively a series expansion for the 
photon number density at the cloud surface: E*(l) — 



£ *[i-tr7(2VO- 



Thus the leading-order approxima- 



tion reduces to the statement that a small cloud blocks no 
radiation in any direction, so E*(l) = Eq, the unatten- 
uated value. The next-order correction accounts for the 
small fraction of photons that are blocked at the cloud 
surface. 

Case 3: Large Clouds. Our final limiting case is that of 
a cloud so large that the transition from atomic to molec- 
ular gas occurs in a thin layer at the cloud surface, so 
that the cloud's curvature is negligible. Before proceed- 
ing we note that this case is not the same as the case of a 
one-dimensional slab subject to a unidirectional beam of 
radiation that we analyzed in § [3] The difference is that 
here the radiation field is isotropic, so it has an angular 
dependence that can vary with depth within the cloud. 
For this reason the large cloud limit with an isotropic 
radiation field is a two-dimensional problem even if the 
cloud is a semi-infinite slab. To analyze this case we per- 
form a series expansion around the limit tr — * 00, but 
with j/h 2 tr finite, so that there is a finite optical depth 
to the molecular region. If ip > 2 then a dust-dominated 
zone exists, and we solve this problem by starting from 
equations (|All|) and (|A12|) and series expanding 770 and 



r/i to first order in t r 



r)o(xd,xn 2 ;Tn) = 
where Td 



Doing so gives 

E 2 (r d ) | e- Td 4 
2 

E 3 (Td) 



TdE 2 (r d ) 



VdTR- 

1 



Equation ([33 
E 2 (r d ) , e" T ' 



4r R 
TdE 3 (r d ) 
2r R ' 
therefore becomes 
+ TdE 2 (r d ) 



4r R 



(52) 
(53) 

(54) 



to first order in , which is straightforward to solve 
numerically to determine Td for a given ip and tr . Alter- 
nately, one may obtain a purely analytic expression by 
dropping the 1/tr correction term. In this case the equa- 
tion becomes E 2 (rd) = 2/ip, which has the approximate 
solution 

r d w 0.83 ln(0.27/> + 0.6); (55) 

this expression is accurate to better than 2% for 2 < ip < 
100. Since Xd — 1 — id/rR, this fixes Xd- Similarly, once 
Td is known it is straghtforward to solve equation (|37p to 
first order in t,7 to obtain 



xn. 



1 



Td 



1P 



E 3 (r d ; 



(56) 



For ip < 2, there is no dust-dominated zone, and 
we must instead solve equations (|43|) and (|44| in the 
limit tr — > 00. We note that for a very large cloud 
E*(i) — E£/2 because the cloud blocks half the sky, and 
it therefore follows immediately from the definition of a e 
(equation [55)1 that 



1~cR 



2tr 

iP ' 



(57) 



i.e. that the effective molecular opacity is a factor of two 
larger than its free-space value because the radiation in- 
tensity at the cloud surface has half its free-space value. 
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Fig. 5. — The plot shows the solution in the large cloud limit. 
The curves shown are the dust optical depth to the fully molecular 
region tjjt = 1r(1 — xn 2 ) (solid curve), the dust optical depth to 
the point of dust-molecular absorption equality = Tr,(1 — 
(dashed curve to the right of tf) = 2), and the ratio of the effective 
molecular opacity to the dust opacity t^/tr (dashed curve to the 
left of i/) = 2). The dotted vertical line at vb = 2 indicates the 
boundary between the presence and absence of a dust-dominated 
zone in the weak radiation limit. 

Similarly, the flux is F*(l) = F */4 because half the sky 
is blocked and the radiation direction is random over 
the other half. Using this boundary condition to inte- 
grate the one-dimensional transfer-dissocation equation 
(|10j) with dd = and the molecular absorption rate mul- 
tiplied by 4> then gives 



xh 2 = 1 - 



4t r ' 



(58) 



We verify that these intuitive arguments in fact give the 
correct leading order terms in the series expansion in 
Appendix O Thus we have the limiting solution to first 
order in 1 for both ip < 2 and ip > 2. We illustrate 
this solution for tr — > oo in Figure O 

4.5. Numerical Solution and Analytic Approximation 

We now proceed with a numerical treatment of the gen- 
eral case. We solve equations (|33| and (|37j) . or (j43|) and 
(j4~4")l , on a grid of points in the (tr, -0)-plane, and plot the 
results in Figure[6] In addition to plotting xh 2 and either 
or t c r, we also show two derived quantities of inter- 
est. The first is the dust optical depth to the molecular 
transition along a radial trajectory, thi = nR(l — xn 2 )o~d- 
We may think of this as the Hi "shielding column" times 
the dust cross section. The second is x# , which is the 
fraction of the cloud's volume that is within the predom- 
inantly molecular region. 

The general behavior of these curves can be understood 
intuitively. If one fixes the cloud density n and dust opac- 
ity (Td, then as the cloud radius increases so does tr, and 
for fixed external radiation field ip the molecular transi- 
tion moves outward, but the Hi column to that transition 
approaches a constant value. Similarly, at fixed cloud 
size and hence tr, increasing the external illumination ip 
raises the amount of atomic hydrogen that is required to 
shield the molecules. Thus ih 2 drops when ip increases 
at fixed tr. 

These curves also enable us to determine under what 
circumstances dust makes a significant contribution to 
shielding the gas, a subject that has been discussed con- 
siderably in the literature (e.g. Ivan Dishoeck fc Black! 



19861 iDraine fc Bertoldil Il996) . To evaluate the impor- 
tance of dust, we can consider how the molecular and 
atomic volumes change as <7d — * 0. In terms of our pa- 
rameters, this amounts to taking the limit as x —> an d 
tr — * 0, but the ratio x/tr remains constant. Graphi- 
cally, this is equivalent to sliding toward the lower right of 
Figured! along a trajectory that is close to a line of slope 
unity - it is not precisely a line of slope unity because of 
the slight non-linearity of the relationship between ip and 
X- We plot the factor by which dust shielding changes 
the radius of the molecular zone or the radial path length 
through the atomic zone, whichever is larger, in Figure 
[7l As the figure shows, dust shielding changes the ra- 
dius of the molecular zone by a factor of ~ 2 only when 
ip is of order unity or larger or when tr is very large, 
which is about what one might expect. Dust shielding 
can affect the size of the molecular region even if there 
is no dust-dominated zone because we have allowed for 
dust absorptions even in the molecular shielding region. 
However, the effect is at most tens of percent. A larger 
change is possible only if the radiation field is intense 

enough to create a dust-dominated zone (i.e. tp>l) or 
the cloud is so large (i.e. tr 3> 1) that even weak dust 
shielding becomes significant because it attenuates the 
radiation exponentially rather than in a powerlaw fash- 
ion as do the molecules. 

We can deduce approximate analytic fitting formulae 
for ih 2 by interpolating between the solutions for the 
limiting cases. The following fitting formulae are reason- 
ably accurate and can be evaluated with no numerical 
iteration: 



■4> 4 + (3r n ) 4 ■ 

max(x s ,x^), 



xu 2 



-3/2 



-3/2 
— 3\ — 1/ ' 



-2/3 



0, 



Ip < 1 
1< 1p < 1p h 

iP h <iP<2 (59) 

2 < -0 < e TR 
e TR < ip 



where ipb = (1-4 + 2tr)/(1.4 + tr) is the approximate 
value of ip at the boundary between the existence and 
non-existence of a zone of dust-dominated opacity, 



JOL. 

\V3Tp 



2V^t| 

5tf>' 2 



x: = < l 



i 

4ri 



1536 
25tr(tr+2) 3 



In ( 1 -t- zT b 

1/2 



5/2 



tp < 1 
1< Ip < lp h 

ifa < i> 

is the approximate value of xr 2 m the strong radiation 
(for ip > 1) or small cloud (for ip < 1) limits (note that 
we only evaluate x s when ip < e TR ), and 



x l 



^max 



1-7^,0 

4t r ' 



1 

4t,:' 



1- T d 



In (1 



2tk 



4+2r d / tr 



Ip <1pb 

ipb < ip < 2 
2 < V 



.( 61 ) 

is the approximate value of ih 2 m the large cloud limit. 
Here 

T d = E^ 1 w 0.83 ln(O.20 + 0.6) (62) 



11 




Fig. 6. — Contours showing the solution as a function of tr and ifi for the structure of the PDR in the two-zone approximation. The 
values shown are, clockwise from the upper left, xh 2 , £d or Tr/t"cR.! thi, and £g 2 - The hatched region is the region in which there is no 
primarily molecular part of the cloud. The dotted line indicates the boundary between the presence and absence of a dust-dominated zone. 
In the panel labelled xa, ra/r e R, the contours above the dotted line indicate the value of Xd, those below it show the value of Tn/r e n_, and 
on the dotted line both of these quantities are exactly 1.0. We caution that the contours for ih 2 = 0.1 and x^ 2 = 0.3, and for ig 5 = 0.1, 
should be regarded as giving upper limits on ih 2 : n °t precise estimates - see the discussion in § 1 1.71 




Fig. 7. — Contours showing the factor by which dust shielding 
changes the atomic or molecular volume at a given (tr , if)) . The 
quantity plotted is max[i H2 /iH 2 ,iicl. (1 — a; H 2 ,nd)/(l — X H 2 )L where 
the subscript "nd" indicates the value with no dust shielding in the 
limit, i.e. in the limit <7d — » 0. The quantity plotted is therefore 
the fractional amount by which dust shielding increases the radius 
of the molecular zone or decreases the radial path length through 
the atomic zone, whichever is larger. 

is the approximate optical depth to the dust-molecular 
opacity crossover in the large cloud limit when dust 
shielding is important, 



is the value of tr at the dust-no dust boundary for a given 



tp, and for convenience we have used the approximation 
Es(x) « e~ x /(2+x). Note that these approximations can 
fail if one is very near the dust-no dust boundary because 
the approximation Tb ~ lA(ip— 1)/(2— ip) is insufficiently 
accurate; in this case one may still use the approximate 
expressions by replacing Tb with a more accurate value of 
tr, on the dust-no dust boundary computed as described 
in § 1431 

Figure [5] shows the error in our analytic approximation 
as a function of tr and ip, where we define the error 
as \xn 2 - x H2 ,approx|/ max(x H2 , 1 - xn 2 ) and xn 2 is the 
solution obtained by numerically solving the appropriate 
equations. As the plot shows the fitting formulae are 
generally good to the ~ 10% level, as good as the two- 
zone approximation itself. The maximum error over the 
range 0.1 < r R < 100 and 0.1 < V < 100 is 24%, and 
occurs near tr = 1.8, ijj — 1.9. 

We can also obtain an even simpler approximation for- 
mula if we specialize to the case where there is no or 
almost no dust-dominated zone and ccr 2 ^ 0.5, which we 
show in Paper II is the most common case in nearby 
galaxies. Consider equation (|4"4"| . which describes the 
surface flux for the case of no dust. The two terms on 
the LHS represent the contributions to the flux from rays 
that do and do not pass through the cloud, respectively. 
If the cloud has a significant molecular core, xh 2 ^ 0.5, 
then only for a small range of angles do rays pass through 
the cloud but not strike the opaque molecular core, and 
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Fig. 8. — Error in the approximate analytic fit given in § 14,51 
as a function of tr and V- The shaded regions indicate errors 
below 2.5% (no shading), 2.5-5%, 5-10%, 10-20%, and > 20% 
(darkest shading). The maximum error is 24%. The dotted lines 
show the boundaries of our different approximation regions: ip < 1, 
1 < ip < Vbi V'b < < 2, and 2 < ip. The hatched region is 
ip > e r R, where there is no predominantly molecular core. Note 
that the error jumps at the X = 1 an d X = 2 lines because the 
fitting formula is slightly discontinuous there. 

thus the first term on the LHS is small in comparison to 
the second. For convenience we define 



5 = 



1 - (1 - 2^H 2 T eR )e 2 ^ T '* 



2r 2 



which enables us to rewrite equation (|44[) as 
3 , 3^ „ « , 3^M 



(64) 



(65) 



where 5 is a small, positive number. Now consider how 
S varies with tr: we show in § 14.41 that for either small 
or large tr, to first order t r oc T^/ip. If we consider the 
series of expansion of 5, this implies that 6 approaches 
a constant at small tr and varies as V ,2 / t r for large tr. 
To generate our approximation we adopt an intermediate 
scaling 



TR 



(66) 



where a is a constant to be chosen to optimize the ap- 
proximation. This gives 



1 - 



3ip 



4(t r + aip) 



(67) 



For the choice a — 0.2, equation (|6~7| agrees with the 
numerical solution to equations (f33|) and (j37|) or (|43|) and 
(|44|) to better than 15% whenever V> < 3 and equation 
(|67[) gives rr^ 2 > 0.15. A corresponding approximate 
formula for the optical depth through the atomic layer is 



4tr - a'tp 



(68) 



with a' 



4a. For a = 0.2 this expression agrees with 



the numerical solution to better than 15% for tp < 3 
whenever equation (|67|) gives i| 2 > 0.1. 

4.6. Comparison to the One- Dimensional Case 

Now that we have solved the spherical case, we are in 
a position to compare to the case of a one-dimensional 
beam of radiation striking an infinite slab that is often 



treated in the literature. This will allow us to deter- 
mine when this approximation yields reasonably accu- 
rate results, and when it is it gives significantly different 
results. Figure [9] shows a comparison between our so- 
lution with isotropic radiation and varying cloud sizes 
versus the most common approximation in the litera- 
ture: an infinite cloud and a beam of radiation whose 
photon number density is half the free-space value. In 
the calculations for finite clouds and isotropic radiation 
we end each curve at the value of ip for which the fully 
molecular region vanishes. For the beamed radiation and 
infinite cloud case, we use the analytic solution described 
in §0 

As the plot shows, when tr <C ip, the one-dimensional 
slab approximation can produce significantly different es- 
timates of the depth of the dust shielding layer than 
does our higher-dimensional approach. The difference 
becomes larger as we consider smaller clouds. For tr > 
1 the slab approximation generally underestimates the 
depth of the atomic layer by tens of percent, primarily be- 
cause it assumes neglects the photodissociation provided 
by non-radial rays. Even for a cloud that is infinitely 
large, tr = oo, this difference between an isotropic radi- 
ation field and a beamed one can be significant at mod- 
erate ip because even though there are no rays reach- 
ing a given position from the "back side" of the cloud, 
pi < 0, when the radiation field is isotropic there are still 
non-radial rays that raise the photodissociation rate at 
a given position above what it would be in the purely 
beamed radiation field of smaller intensity. 

For tr < 1 the sign of the error depends on ip. When 
the radiation field is weak, the slab approximation also 
underestimates the depth of the atomic layer, for the 
same reason as when tr > 1. When the radiation field 
becomes strong, though, the sign of the error reverses, 
although as we discus in § 14.71 our fiducial model is of 
limited accuracy for small tr and large ip. 

Physically, clouds of a wide range of sizes and densities 
are of course present in the ISM. For the atomic envelopes 
of GMCs in the Milky Way, a typical density is n ~ 30 
cm~ 3 and a typical dust cross-section is 10 -21 cm 2 , so 
that ~ 10 pc of path provides an optical depth of about 
1. Since these envelopes are a few tens of pc in size, a 
typical one might have tr of a few, in which case the slab 
treatment underestimates the true size of the envelope at 
the tens of percent level. In low-metallicity galaxies with 
low molecular fractions, however, the error is likely to be 
much worse because tr will be significantly smaller. 

4.7. Uncertainties in Spherical Geometry 

We have shown that in the case of a one-dimensional 
beam of radiation impinging on a slab, the two-zone ap- 
proximation is capable of determining the neutral hy- 
drogen shielding column to better than ~ 50% accuracy. 
This characterizes the level of error imposed by most 
of our physical assumptions. However, in spherical ge- 
ometry we have an additional uncertainty, imposed by 
the fact that we must assign an effective optical depth 
to rays that pass at arbitrary angles through the region 
where molecular shielding dominates, but the gas is not 
yet fully molecular. In particular, rays from the "back 
side" of our cloud, those with fj, < 0, contribute to the 
energy density and flux throughout the cloud. Such rays 
are absent in the case of an infinite planar cloud, because 
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TRX d 

3ip 



Fig. 9. — The upper panel shows the dust optical depth to 
the point where the gas becomes predominantly molecular tjji for 
an isotropic radiation field of normalized intensity tjj and various 
values of tr (solid lines) and for a unidirectional radiation field of 
normalized intensity ip/2 (dashed line). The lower panel shows the 
fractional difference between the results for finite tr, and isotropic 
radiation and for an infinite slab illuminated by unidirectional ra- 
diation, defined as Difference = [thi iso — T HI beam] /"Til, iso- In all 
cases the curves for finite tr end at the value of tp for which fully 
molecular region disappears (-0 = e TR ). 



rays with /i < pass through the fully molecular region 
and are therefore infinitely attenuated. 

In this section, we seek to determine how much addi- 
tional uncertainty is introduced into our calculations in 
spherical geometry by this complication. To do so, we 
note that in our treatment above we assume that the 
opacity in the molecular shielding region will be roughly 
equal to that at its surface, i.e. that it does not rise 
sharply until one is very close to the transition to fully 
molecular gas. The represents a minimum attenuation 
along [i < rays. To check the importance of that as- 
sumption, we consider an extreme assumption in the op- 
posite direction: that all rays with /i < are infinitely 
attenuated. This assumption is obviously unphysical, 
since if it were true then the transition to fully molec- 
ular would occur as soon as self-shielding began to dom- 
inate over dust shielding. However, it provides a worst 
case with which we can compare our fiducial model as a 
way of characterizing our uncertainty. Since, as we shall 
see, the value of xn 2 that we obtain by making this as- 
sumption is always smaller than what we obtain for the 
fiducial case, we may regard the fiducial case as giving 
upper limit on xh 2 and the case we calculate here as 
giving a lower limit. 

If we take the intensity along rays with /i < to be 
zero, this is equivalent to replacing r]o(xd, ^h 2 i T R.) with 
770 (xd, Xd] tr) in the equations derived in § 14. II and § 14.2) 
and similarly for rji- Doing this and simplifying gives 



1 



xr 2 

Xd 



(70) 



1 



(71) 



r] (xd,x d ;T R ) 



1 



(69) 



for tp > 2 (dust shielding is significant) and 

3tP_\ 1/3 
4r R J 

for -0 < 2 (dust shielding is not significant). Note that in 
this case when dust shielding is negligible it is possible to 
solve the equations analytically, which we have done to 
obtain equation (|7ip . It is immediately obvious from this 
equation that the fully molecular region vanishes in the 
region 4tr/3 < tp < 2. We plot the solutions to equations 
(|6"9"|) - (f7"Tj) in the (tr, tp) plane in Figure [TOl and we show 
the difference between this solution and our fiducial one 
in Figure [TTl 

As the plots show, the difference between the two mod- 
els is negligibly small over most of parameter space; 
there is a significant difference only in the region roughly 
bounded by the curves tp < e TR , tp > 4tr/3, and tp < 2. 
Alternately, we can phrase these constraints in terms of 
values of . This is particularly useful for tp < 1 , where 
contours of both constant xh 2 and constant uncertainty 
are both straight lines corresponding to fixed V'/tr. The 
10%, 50%, and 100% uncertainty contours in x^ 2 , shown 
in the right panel of Figure [11] correspond to values of 
= 0.47, 0.28, and 0.20, respectively, as computed us- 
ing our fiducial model and shown in Figure [TTl Since we 
have already established that the two-zone approxima- 
tion is uncertain at the tens of percent level, the geomet- 
ric uncertainty is probably only dominant when tp < 1 
and our predicted molecular volume fraction is less 
than about a quarter. If tp > 1, the errors at a given 
value of Xft 2 are considerably smaller, so the geometric 
uncertainty is not important except for clouds with very 
small molecular fractions. At such molecular fractions, 
one should interpret our fiducial case as giving only an 
upper limit on the molecular content of a cloud. 

The significant geometric uncertainty for such clouds 
is not surprising, since these clouds are near the limit 
of having no molecules at all. For them any change in 
our physical assumptions that increases or reduces the 
amount of shielding even a small amount produces a sig- 
nificant change in the results. Indeed, our fiducial cal- 
culation shows some unphysical behavior in this regime, 
in that we find that for tp < 1 there is no finite value of 
tr for which the fully molecular region vanishes and the 
cloud remains atomic throughout. This seems unlikely, 
given that for a chosen value of tp and very large tr the 
thickness of the atomic region approaches a finite value; 
one would expect that clouds much smaller than this 
should be atomic throughout regardless of their shape, 
and indeed under the assumption of infinite attenuation 
for backside rays that we make in this section, the molec- 
ular core always vanishes at some finite value of tr for 
any finite tp. 

Fortunately, as we discuss in § 14.81 for realistic pa- 
rameters describing giant atomic-molecular complexes, 
we generally have tp > 1, tr > 1, and, as we show in 
? 14.81 x^ 2 > 0.5, and in this part of parameter space the 
uncertainty introduced by the /i < rays in spherical 
geometry is < 10%. 

4.8. Example Calculations 
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F ig. 10. — Same is Figure[6] but assuming infinite attenuation along rays that enter the molecular absorption-dominated region (equations 
[MED- 
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Fig. 11. — The plots show the difference between thi and computed under the fiducial assumption that the opacity throughout 
the molecule-dominated region equal to that at its surface (§ 14, H and § 14.21 1 and under the assumption of infinite attenuation in this region 
(§ 14. 71 . The difference is defined as [thi, fiducial — T Hl,attcnuatcdl/ T Hl,fiduciali an( i similarly for ^j^- The hatched region is the region in 
which there is no predominantly molecular core under either assumption; the difference in this region is obviously zero. 



Here we provide examples that illustrate the use of 
our analytic approximations for PDR structure. Since 
these calculations are intended to be illustrative rather 
than to analyze real situations (which we will discuss in 
the next paper in this series), we choose parameters to 
yield examples that span the possible combinations of 
parameters without worrying how well they agree with 
observations. As a first case, consider a "typical" Milky 
Way cloud, with ad = 1.1 X 10~ 21 cm 2 , n — 30 cm~ 3 , 



K = 3x 10" 17 



cm s 



and cE* = 10 8 cm' 



L . Note 



that we have used a value of Eq somewhat larger than the 



solar neighborhood value because most molecular clouds 
are closer to the galactic center, where the radiation field 
is more intense. This combination of parameters gives 
X = UssCT d cE^/(nn) = 12.2 and V = x(2.5 + X )/(2.5 + 
%e) = 6.3. Now consider a giant molecular cloud complex 
with a radius R = 50 pc, which gives tr = na^R = 4.64. 
Since ip > 2, this cloud has a significant dust shielding 
zone, and since ip < e TR it also has a fully molecular core, 
as we expect. Using the approximation equations (|59l) 
(jnj) for this case, we find Td = 0.52, xg — 0.74, x s = 
14.1, and an approximate value of ih 2 ~ 0.74. (Note 
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that although xn 2 is strictly less than unity, it is possible 
for x s to be larger than unity because we only retained 
a finite number of terms in the series expansion used 
to generate it. However, due to the way x s and xi are 
combined, our approximate expression for xh 2 is always 
less than unity.) Numerical solution for these parameters 
gives ih 2 = 0.70. Such a cloud is 34% molecular by 
volume, and is shielded by an atomic column that is 15 
pc deep and has a column density of Nm — 1.4 x 10 21 
cm" 2 from the edge the cloud to the edge of the molecular 
zone. 

Now consider moving this cloud to a point farther out 
in the Galaxy where the ambient FUV radiation field is 
weaker, so that all cloud parameters remain the same 
but now cEq — 2 x 10 7 cm -2 s _1 , a factor of 5 below our 
previous value. In this case we have the same tr, but 
X = 2.44 and -0 = 1-43. For this r R we have tpb = 1-40, 
slightly smaller than i/j, so this cloud just barely still has 
a dust-dominated zone. Evaluating equations (|59 p - ([ST|) . 
we have X£ = 0.87 and a; s = 95, so get an approximate 
value xu 2 ~ 0.87; the numerical solution is ih 2 = 0.92. 
Thus moving the cloud to this reduced-radiation envi- 
ronment raises the molecular volume fraction to 77%, 
and reduces the Hi shielding column to a layer 4 pc deep 
containing a column of Nm = 3.7 x 10 20 cm" 2 hydrogen 
atoms. If the cloud were slightly denser, n = 40 cm" 3 
instead of 30 cm' 3 , then tr would increase to 6.18 from 
4.64, and ip would decrease to 1.11 from 1.43. Since 
ipb = 1.82 in this case, the cloud would be dominated 
by molecular absorption throughout. Evaluating the ap- 
proximation equations gives xi = 0.95, x s = 0.35, and 
xn 2 « 0.95; the numerical solution is xh 2 = 0.95. Thus, 
the increase in density would slightly increase the molec- 
ular volume and the column density through the shield- 
ing layer to 87% and Nm — 5.2 x 10 20 cm" 2 , respectively. 

Finally consider a cloud in a low-pressure dwarf galaxy 
with a very low star formation rate, so the cloud has 
lower density and metallicity than a Milky Way cloud, 
n=W cm" 3 , cr d = 2.2 x 10" 22 , and K = 6x 10" 18 , and 
is exposed to a lower level of radiation, cEq = 10 6 cm" 3 
s _1 . We keep the cloud radius unchanged. This cloud 
has tr = 0.31 and tp — 0.29, which from our approximate 
formulae gives xi = 0.77, x s = 1.4, and xn 2 rs 0.77. The 
numerical solution is xh 2 — 0.73. This cloud would be 
38% molecular by volume, and would have a shielding 
layer of Nm = 4.2 x 10 20 cm -2 , 14 pc deep. 

5. SUMMARY AND CONCLUSION 

In this paper we develop an approximate analytic so- 
lution to the problem of determining the size of the PDR 
that bounds a cloud of gas embedded in a dissociating 
background radiation field. This is a reasonable approx- 
imation to the problem of finding the location of the 
transition between the atomic envelope and the molecu- 
lar core in a giant atomic-molecular cloud complex, such 
as those which contain the bulk of the molecular gas in 
the Milky Way. 

We show that the location of the transition is deter- 
mined by two dimcnsionless parameters. These are tr, 
the dust optical depth through the cloud, and x, the ra- 
tio of the rate at which dissociating photons are absorbed 
by dust grains to the rate at which they are absorbed by 
H2 molecules in the absence of any shielding. We may in- 
tuitively think of these parameters as characterizing the 



size of the cloud and the intensity of the radiation field 
to which it is subjected. Within this parameter space we 
identify two critical curves, which define the boundaries 
at which a fully molecular region in the cloud center ap- 
pears, and at which dust shielding begins to contribute 
significantly to the shielding of H2 molecules. We develop 
the equations that determine the sizes of the molecular 
and atomic regions in this parameter space, and we pro- 
vide an approximate analytic solution for them (equa- 
tions [SH] - |5T] and equations [57] and |6"5]) . Our solutions 

are accurate to tens of percent for clouds that are > 20% 
molecular by volume, and provide upper limits on the 
molecular content at this accuracy for clouds with lower 
molecular content. Using this formalism we find that for 
typical giant atomic-molecular complexes in the Milky 
Way x ~ 1) which indicates that dust shielding and self- 
shielding each make order unity contributions to deter- 
mining the location of the atomic-molecular transition. 

Our work shows that the procedure of determing the 
structure of PDRs by treating them as semi-infinite slabs 
illuminated by unidirectional beams of dissociating radi- 
ation is a reasonable approximation for extremely opaque 
clouds, but that it fails badly for small clouds or weak 
radiation fields, i.e. in cases where the transition from 
atomic to molecular is sufficiently far into the cloud that 
the cloud's curvature cannot reasonably be neglected. In 
such cases the slab approximation can either overesti- 
mate or underestimate the size of the atomic layer by 
factors of order unity, depending on the particular pa- 
rameters of the cloud and the ambient radiation field. 

The development of an analytic model for the structure 
of the atomic envelopes of finite molecular clouds opens 
up the possibility of developing a more general theory 
of the atomic to molecular ratio in galaxies. In a galaxy, 
the mean interstellar radiation field and the conditions in 
the atomic portion of the atomic interstellar medium are 
determined by the star formation rate, which determines 
the abundance of young, hot stars. In turn, the star for- 
mation rate depends on the fraction of the ISM of that 
galaxy in molecular form, and therefore available for star 
formation. At some level, therefore, star formation in 
galaxies much be a self-regulating process, with the for- 
mation and dissociation of molecular clouds representing 
one step in that regulation. Developing a simple model 
for how the molecular fraction in a cloud is determined 
by its properties and those of the ambient radiation field 
represents a step toward a complete theory of the star 
formation rate. In future work, we plan to develop this 
theory further by applying the model demonstrated here 
to molecular clouds in galaxies. 
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APPENDIX 
A. EVALUATION OF no AND m 



Here we evaluate the two functions 



1 f 1 

Vo{x d , xh 2 ] tr) = - / d/x exp(-r R £) 

1 f 1 

m(xd,xu 2 ;T R ) = - / d^/icxp(-r R ^). 



(Al) 
(A2) 



To evaluate the integrals, we change the variable of integration from [i to £. Using the definition of £ (equation 
we find that 

(i - Xd )(i + Xd ) - e 



and 

Making the change of variable, we find that 

•l—xa 



dfi {1 - x d )(l + x d ) + £ 2 



(A3) 
(A4) 



l — x d 



m(xd,XH 2 ;Tn) - ^ , 



l 

8x1 



d£ exp(-£r R ) + (1 - x d )(l + a*) / d£C 2 ex P (-£r R ) 



1— x d rl—Xd 

^eexp(-^r R ) - (1 - x d ) 2 (l + x d ) 2 / ^r 3 exp(-^r R ) 

?h 2 ^?h 2 



where 



& 2 = £(z d ,MH 2 ) = \ /L-xl + xlfil -x d /x H2 - 



(A5) 



(A6) 



(A7) 



Of these integrals, the first one for 770 may be evaluated directly, while the first one for r/i may be evaluated by parts. 
The second integral on each line may be evaluated via the identity 



dxx- n e~ ax = dxx- n e~ ax - dxx~ n e' ax 

X J Xq J X! 



■ X, 



.1-n 



duu e — Xi / aw e 
1 Ji 



--x n E n (ax )-x 1 n E n (axi), 



(A8) 

(A9) 
(A10) 



where in the second step we made the change of variables u — x/xq and v — x/xx, and E„ is the exponential integral 
function of order n, defined by E„(x) = J, t~ n e~ xt dt. Using this identity gives 



7]o(Xd,XH 2 ;T R ) = 



4i d 



TR 



+ (l+a?d) 



E 2 (2MTr) - — ^E 2 (£h 2 t r ) 



77i(x d ,XH 2 ;rRj 



1 j (l + CH 2 r R )e-« H 2-R - (1 + 



8^1 



ydTRje 



£h 2 



+ (1 + x d ) 



E3(ydTR)-^-E 3 (?H 2 r R : 

?H 2 



(All) 
(A12) 



where y d = 1 — x d . Note that exponential integrals obey the recurrence relation nE n+ i(i) = e x — xE n (x), so we 
could alternately have written these in terms of Ei(a;) or the classical exponential integral Ei(x) = — Ei(— x). 

B. SOLUTION BY SERIES EXPANSION IN THE STRONG RADIATION LIMIT 



Here we solve equations (|33j) and (I37|) in the strong radiation limit, i.e. xh 2 *C 1, x d <C 1, and /in 2 + 1 *C 1, by 
means of series expansion. Let (3 = x d T R and 7=1 + /iu 2 • Then we have 



1 + M/3 + ( ± + — ^— — M 2 ^) /? 2 + (V^ + ~7 — "M 3 ) /3 3 + 0(/3 4 ) 



2r R 



2r R 



6r R 



(Bl) 
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Note that we have retained terms out to order /3 3 . We shall see below that this is required for a consistent solution. 
Using this expansion in the integrals 



1 f 1 
1 f 1 

rn(xd,XH 2 ;T R ) = - J c^e m< = 



and expanding in powers of 7, we find 
?7o(^d,^H 2 ;rR) = — 



77i(x d ,XH 2 ;T R ): 



2 - 7 + 7/5 + + 0(/3 4 ) + 0( 7 /3 2 ) + 0( 7 2 ) 

OTr 

2^ „ r R + 2 3 /->/|C}4\ 1 r\l a2\ , .n/ 2\ 



3/3 + 7-7/3+' 15 



-^ + 0(/3 4 ) + 0( 7 /^) + 0( 7 ') 



If we similarly expand the right-hand sides of (|33[) and (|3T[) in powers of (3 and 7, we obtain the two equations 

= i + 0(/3 4 ) + 0( 7 /3 2 ) + 0(7 2 ) 

:A + 0(/? 4 ) + 0( 7 /? 2 ) + 0(7 2 )- 
3^ 



2-7 + 7/3+^/3 2 
3r R 

2 - . - ra + 2 „ 3 



3/3+7-7/3+ 15 



(B2) 
(B3) 

(B4) 
(B5) 

(B6) 
(B7) 



If we combine these two equations by eliminating the common factor e rR / - 0, we obtain an equation for the relationship 
between and 7: 



7 7/3 r R + 2 2 = 3 /7 
2 2 + 6r R P 2 V/3 



- ^7 + + (/3 3 ) + W) + O 



(B8) 



The only way for this equation to have a consistent solution in which the orders on both sides balance is if 7 is of 
order /3 3 . In this case the leading order on both sides is j3 2 (an order we retained only be performing the expansion in 
equation IB II to order /3 3 ), and balancing the leading order terms gives 



7 



2tr + 4 3 
45r R 



Since we now know the order of all terms, we can solve equation (|B6[) to leading order to obtain 



Xd 



Similarly, we know that 



£ 

TR 



XR 2 



G 



tr(tr + 2) V tp 



1/2 



-2L) =1-M H2 =27 + 0(7^ 

Xd 



Substituting (|B9j) for 7 and (IBlOjl for a;d and re-arranging, we obtain to leading order 

1536 



xh 2 



25tr(tr + 2) s 



1 / e r R \5, ! 



(B9) 

(BIO) 
(Bll) 

(B12) 



C. SOLUTION BY SERIES EXPANSION IN THE LARGE AND SMALL CLOUD LIMITS 



Here we solve equations (|4"3")l and (|4"4"]l by series expansion in the limits r R — > and r R — > 00. We approach this 
problem by defining — t c r/tr, so that with some rearrangement the equations are 



1 _ e 2/3MH 2 r R + 2/3tr _ _ Tr = Q 



(1 - 2f3 m ^)e^^ _ 1 + ^2 _ ^3 ^ _ (1 _ ^3/2 



= 0. 



For the case tr — > we then let 



/3 = /3 +/3lTR + #2T R + 



Mh 2 = Mo + Mi r R + M2T R H . 

Expanding equations (IC1[) and ()C2[) to leading order in r R and re-arranging gives 

2^(1 -Mo) -7 = 

2/3 2 (l- M 2 )=0, 



(CI) 
(C2) 

(C3) 
(C4) 

(C5) 
(C6) 
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which has the solution (Xa = — 1, /3q = l/ip. Using these values and continuing the expansion to the next order, we 
obtain 



2/3l + i±^ =0 



4/xi 

2 



so Hi = and (3i = l/(2-0 2 ). Continuing to one more order, we have 



(C7) 
(C8) 

(C9) 
(CIO) 



-1 + 6^=0, 

so ju,2 — l/(6"0 2 ) and /3% — l/(4-0 3 ). It improves the accuracy of the approximation for xn 2 significantly at small ip to 
include one more order, so we do so: 



4 ^3 + ^T4--M3=Q 



1 2 

4 8 

ip 2 ^ 5ip 5 



= 0, 



so = 2/(5V> 3 ) and f3 3 = 7/(6O0 4 ). Therefore to order r R we have 

T e R 1 



Xh 2 



1 




1 — + 


7r R 


* + 


2^ 2 " 


4^ 3 




-1 - 




, 2r| 
















f 2\/3 

5-0 




-/4 2 


V30 



For the case r R — ► oo, we let 



/3 = /3o + /3-i/ 2 t r 1/2 +/3-it r 1 + .-. 

i"h 2 = Vo + M-i/27-r 1/2 + Z-'-iTr 1 H , 

and if we expand equations (|C1[) and (|C2[) in powers of t^ 1 then the leading order equations are 

4 

2/3 - T = 

■0 



3V> 



i - (i - 14)*'* 



= 0. 



Therefore (3q = 2/ip and /Xq = 0. Continuing the expansion to the next order, 



2/3- 



1/2 : 



16 

0' 



2 







= 0, 



so 0-i/2 = and ^_i/ 2 



'■0/2. Continuing one more order, 
l + 2/3_i=0 
16^2 



0> 5 / 2 



M-i=0, 



so /3_i 



T/2 and = Thus to order t^ 1 in the limit r R — > oo we have 



T~cR 


2 




TR 






m 2 




/ 1/> 

2r R 


IH 2 


= 1 - 


4r R ' 



(Cll) 
(C12) 

(C13) 
(C14) 
(C15) 

(C16) 
(C17) 

(C18) 
(C19) 

(C20) 
(C21) 

(C22) 
(C23) 

(C24) 

(C25) 
(C26) 
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